In [2]:
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from libpysal.weights import Queen
import esda
from esda.moran import Moran, Moran_Local
import contextily as ctx
from splot.esda import moran_scatterplot, plot_moran, lisa_cluster,plot_moran_simulation
import libpysal as lps
import matplotlib.pyplot as plt
import plotly.express as px
from shapely import wkt
from shapely.wkt import loads
import branca.colormap as cm
import folium
import mapclassify
df = pd.read_csv("merged_MAR16.csv")
gdf = gpd.GeoDataFrame(df, geometry=df['wkt'].apply(loads))
gdf.set_crs(epsg=3857, inplace=True)
gdf = gdf.dropna()
invalid_geometries = gdf[~gdf.is_valid]
print(f"Found {len(invalid_geometries)} invalid geometries")
gdf['geometry'] = gdf.geometry.buffer(0)
gdf['snap_rate'] = gdf['snap_rate'].astype(float)
gdf['chd_pct'] = gdf['chd_pct'].astype(float)
gdf['lowaccess_pct'] = gdf['lowaccess_pct'].astype(float)
gdf['lowaccess_li_pct'] = gdf['lowaccess_li_pct'].astype(float)
gdf['snap_rate'] = gdf['snap_rate'].astype(float)
gdf['total_pop'] = gdf['total_pop'].astype(int)
gdf['centroid'] = gdf.geometry.centroid
centroids_geo = gdf['centroid'].to_crs(epsg=4326)
gdf['longitude'] = centroids_geo.x
gdf['latitude'] = centroids_geo.y
gdf = gdf[["geometry", "snap_rate", "chd_pct", "lowaccess_pct", "lowaccess_li_pct", "geoid20", "SPA_NAME", "pop_ages_65_older_pct", "latino_pct", "unemployed_pct"]]
Found 186 invalid geometries
In [3]:
spa_colors = {
'West': 'red',
'San Fernando': 'blue',
'South': 'green',
'San Gabriel': 'purple',
'Antelope Valley': 'orange',
'South Bay': 'yellow',
'Metro': 'brown',
'East': 'gray'
}
index_colors = {
0: '#ffffcc', # light yellow
1: '#ffeda0', # light orange
2: '#fed976', # slightly darker orange
3: '#feb24c', # orange
4: '#fd8d3c', # dark orange
5: '#fc4e2a', # reddish orange
6: '#e31a1c', # red
7: '#bd0026', # dark red
8: '#800026', # burgundy
9: '#400026', # dark burgundy
10: '#000000' # black
}
spa_colors
gdf = gpd.GeoDataFrame(gdf, geometry=gdf["geometry"])
gdf.set_crs(epsg=3857, inplace=True)
Out[3]:
| geometry | snap_rate | chd_pct | lowaccess_pct | lowaccess_li_pct | geoid20 | SPA_NAME | pop_ages_65_older_pct | latino_pct | unemployed_pct | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | POLYGON ((-118.298 34.263, -118.297 34.263, -1... | 0.000000 | 5.192139 | 0.024308 | 0.007821 | 101110 | San Fernando | 20.926756 | 26.955656 | 7.790493 |
| 1 | POLYGON ((-118.277 34.260, -118.278 34.260, -1... | 2.936857 | 4.704301 | 0.465066 | 0.037391 | 101122 | San Fernando | 21.829971 | 8.693563 | 9.590236 |
| 2 | POLYGON ((-118.278 34.256, -118.278 34.253, -1... | 2.701243 | 5.896306 | 0.000000 | 0.000000 | 101220 | San Fernando | 16.345877 | 42.028152 | 12.653688 |
| 3 | POLYGON ((-118.287 34.256, -118.287 34.253, -1... | 2.894356 | 5.398263 | 0.000000 | 0.000000 | 101221 | San Fernando | 19.515442 | 32.055378 | 4.291146 |
| 4 | POLYGON ((-118.286 34.256, -118.286 34.255, -1... | 6.915629 | 5.410280 | 0.000000 | 0.000000 | 101222 | San Fernando | 16.346153 | 46.118233 | 11.773151 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2480 | POLYGON ((-118.510 34.187, -118.508 34.187, -1... | 44.843048 | 4.046243 | 0.000000 | 0.000000 | 980024 | San Fernando | 16.788321 | 15.328467 | 0.000000 |
| 2481 | POLYGON ((-118.248 33.847, -118.247 33.847, -1... | 6.180470 | 4.300235 | 0.503704 | 0.152189 | 980025 | South Bay | 26.277372 | 4.866180 | 13.592233 |
| 2484 | POLYGON ((-118.281 33.767, -118.281 33.768, -1... | 0.000000 | 4.518828 | 0.000000 | 0.000000 | 980031 | Unknown | 11.006289 | 34.591194 | 0.000000 |
| 2487 | POLYGON ((-117.972 34.037, -117.972 34.037, -1... | 21.279613 | 4.543835 | 0.032253 | 0.014417 | 980035 | San Gabriel | 10.917722 | 77.215187 | 10.263929 |
| 2488 | POLYGON ((-118.079 34.043, -118.079 34.043, -1... | 0.000000 | 5.430255 | 0.128709 | 0.046114 | 980036 | Unknown | 13.530928 | 47.680412 | 5.956113 |
2467 rows × 10 columns
In [4]:
classifier = mapclassify.NaturalBreaks.make(k=5)
# Age score
gdf['age_pct_score'] = gdf[['pop_ages_65_older_pct']].apply(classifier)
gdf[['pop_ages_65_older_pct', 'age_pct_score']].head()
# Hispanic score
gdf['hisp_pct_score'] = gdf[['latino_pct']].apply(classifier)
gdf[['latino_pct', 'hisp_pct_score']].head()
# Unemployment score
gdf['emp_pct_score'] = gdf[['unemployed_pct']].apply(classifier)
gdf[['unemployed_pct', 'emp_pct_score']].head()
# Food Index: Age + Hispanic + Unemployment scores
gdf['priority_index'] = gdf['age_pct_score'] + gdf['hisp_pct_score'] + gdf['emp_pct_score']
gdf['priority_index'].unique()
Out[4]:
array([ 6, 5, 7, 4, 3, 8, 2, 1, 0, 9, 10])
In [7]:
chd_pct_min = gdf['chd_pct'].min()
chd_pct_max = gdf['chd_pct'].max()
snap_rate_min = gdf['snap_rate'].min()
snap_rate_max = gdf['snap_rate'].max()
lowaccess_li_pct_min = 0
lowaccess_li_pct_max = 0.05
lowaccess_pct_min = 0
lowaccess_pct_max = 0.05
priority_min = 0
priority_max = 10
def style_function_snap_rate(feature):
return {
'fillColor': colormap_snap_rate(feature['properties']['snap_rate']),
'color': 'black',
'weight': 1,
'fillOpacity': 0.7
}
def style_function_chd_pct(feature):
return {
'fillColor': colormap_chd_pct(feature['properties']['chd_pct']),
'color': 'black',
'weight': 1,
'fillOpacity': 0.7
}
def style_function_lowaccess_li_pct(feature):
return {
'fillColor': colormap_chd_pct(feature['properties']['lowaccess_li_pct']),
'color': 'black',
'weight': 1,
'fillOpacity': 0.7
}
def style_function_lowaccess_pct(feature):
return {
'fillColor': colormap_chd_pct(feature['properties']['lowaccess_pct']),
'color': 'black',
'weight': 1,
'fillOpacity': 0.7
}
def style_function_spa(feature):
return {
'color': spa_colors.get(feature['properties']['SPA_NAME'], 'white'), # Boundary color
'weight': 1,
'fillOpacity': 0,
"dashArray": "5, 5"
}
def style_function_priority(feature):
return {
'fillColor': index_colors.get(feature['properties']['priority_index'], 'white'),
'weight': 0,
'fillOpacity': 0.7
}
colormap_snap_rate = cm.linear.YlOrRd_09.scale(snap_rate_min, snap_rate_max).to_step(n=5)
colormap_chd_pct = cm.linear.Blues_09.scale(chd_pct_min, chd_pct_max).to_step(n=5)
colormap_lowaccess_li_pct =cm.linear.YlOrRd_09.scale(lowaccess_li_pct_min, lowaccess_li_pct_max).to_step(n=5)
colormap_lowaccess_pct = cm.linear.Blues_09.scale(lowaccess_pct_min, lowaccess_pct_max).to_step(n=5)
colormap_snap_rate.caption = "SNAP Rate"
colormap_chd_pct.caption = "CHD Percentage"
colormap_lowaccess_li_pct.caption = "Low Access Low Income Percentage"
colormap_lowaccess_pct.caption = "Low Access Percentage"
bounds = gdf.total_bounds
center_lat = (bounds[1] + bounds[3]) / 2
center_lon = (bounds[0] + bounds[2]) / 2
m = folium.Map(location=[center_lat, center_lon], zoom_start=8)
geojson_data = gdf.to_json()
#snap_rate_layer = folium.GeoJson(
# geojson_data,
# style_function=style_function_snap_rate,
# name='SNAP Rate',
# popup=folium.GeoJsonPopup(fields=['snap_rate', 'chd_pct', 'lowaccess_li_pct', 'lowaccess_pct'],
# aliases=['SNAP Rate:', 'CHD %:', 'Low Income Low Access %', 'Low Access %'],
# labels=True)
#).add_to(m)
chd_pct_layer = folium.GeoJson(
geojson_data,
style_function=style_function_chd_pct,
name='CHD Percentage',
popup=folium.GeoJsonPopup(fields=['snap_rate', 'chd_pct', 'lowaccess_li_pct', 'lowaccess_pct'],
aliases=['SNAP Rate:', 'CHD %:', 'Low Income Low Access %', 'Low Access %'],
labels=True)
).add_to(m)
#lowaccess_li_pct_layer = folium.GeoJson(
# geojson_data,
# style_function=style_function_lowaccess_li_pct,
# name='Low Income Low Access Percentage',
# popup=folium.GeoJsonPopup(fields=['snap_rate', 'chd_pct', 'lowaccess_li_pct', 'lowaccess_pct'],
# aliases=['SNAP Rate:', 'CHD %:', 'Low Income Low Access %', 'Low Access %'],
# labels=True)
#).add_to(m)
#lowaccess_pct_layer = folium.GeoJson(
# geojson_data,
# style_function=style_function_lowaccess_pct,
# name='Low Access Percentage',
# popup=folium.GeoJsonPopup(fields=['snap_rate', 'chd_pct', 'lowaccess_li_pct', 'lowaccess_pct'],
# aliases=['SNAP Rate:', 'CHD %:', 'Low Income Low Access %', 'Low Access %'],
# labels=True)
#).add_to(m)
spa_layer = folium.GeoJson(
geojson_data,
style_function=style_function_spa,
name='SPA Overlay',
popup=folium.GeoJsonPopup(fields=['SPA_NAME'], aliases=['Service Planning Area:'], labels=True)
).add_to(m)
index_layer = folium.GeoJson(
geojson_data,
style_function=style_function_priority,
name='FI Index Overlay',
popup=folium.GeoJsonPopup(fields=['priority_index'], aliases=['Priority Index:'], labels=True)
).add_to(m)
#colormap_snap_rate.add_to(m)
colormap_chd_pct.add_to(m)
#colormap_lowaccess_li_pct.add_to(m)
#colormap_lowaccess_pct.add_to(m)
folium.LayerControl().add_to(m)
m
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook